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cy^ ■ Abstract 
I. 

We propose a model that leads to the formation of non-equilibrium finite-size domains 
^1 in a biological membrane. Our model considers the active conformational change of the 

inclusions and the coupling between inclusion density and membrane curvature. Two special 
cases with different interactions are studied by Monte Carlo simulations. In case (i) exited 
state inclusions prefer to aggregate. In case (ii) ground state inclusions prefer to aggregate. 
When the inclusion density is not coupled to the local membrane curvature, in case (i) the 
I typical length scale (\/M) of the inclusion clusters shows weak dependence on the excitation 
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■ rate (Kon) of the inclusions for a wide range of Kon but increases fast when Kon becomes 
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. sufficiently large; in case (ii) \J M ^ Kon~^ for a wide range of -fCon- When the inclusion 
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density is coupled to the local membrane curvature, the curvature coupling provides the 
upper limit of the inclusion clusters. In case (i) (case (ii)), the formation of the inclusions is 
suppressed when ^o// (-f^on) is sufhciently large such that the ground state (excited state) 
inclusions do not have sufficient time to aggregate. We also find that the mobility of an 
inclusion in the membrane depends on inclusion-curvature coupling. Our study suggests 
possible mechanisms that produce finite-size domains in biological membranes. 
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Chapter 1 
Introduction 



As shown in Figure ITTH a biological membrane is a self- assembled bilayer composed of 
various types of macromolecules such as lipids, proteins, carbohydrates, cholesterol, 
and other materials. |lj A membrane is thought to be similar to a two-dimensional 
fluid, therefore the aforementioned molecules can diffuse on it. Many experimental 
and theoretical studies have shown that the distribution of these molecules in a mem- 
brane is nonuniform. j2] Instead, some molecules aggregate to form domains with 
typical length scale ranging from tens of nanometers to about one micrometer. j3] 
An important subject of current biomembrane research is to find the mechanisms 
that is responsible for the formation of these finite-size domains. From equilibrium 
statistical mechanism, we know that if there are only short-range interactions be- 
tween the molecules, equilibrium macrophase separation can be induced in a system, 
but equilibrium microphase separation can not be induced if there are no long-range 
interactions in a system. However, up to now there has been no clear evidence for 
such long-range interactions in the membrane. Several possible mechanisms for the 
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Figure 1.1: The configuration of membrane is shown in this figure. 



formation of finite-size domains in biological membranes have been proposed, includ- 
ing (i) the kinetics of spontaneous membrane domain assembly depends on the effect 
of membrane recycling ubiquitous in the living cells |3J E] , (ii) cholesterol can induce 
separate domains of two mixed fatty acids P, . Although the true mechanisms of 
raft formation has not been verified yet, it is possible that different mechanisms are 
involved in different situations. 

Since the activities of proteins in the membranes are important in many physi- 
ological processes including ion transport, signal transduction, et al., in this thesis 
we propose that the activities of proteins can also lead to the formation of finite-size 
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domains in the membranes. To concentrate on the basic mechanism behind activity- 
induced finite-size domains, in our model the system contains only one type of lipids 
an one type of active proteins. These proteins are called inclusions in this thesis. 
The inclusions can change their internal states, and different internal states of in- 
clusions have different conformations and different interaction with other molecules. 
Therefore this work is a natural extension of Ref 0, in which the fiuctuation spec- 
trum and the stability of the homogeneous sate in such systems were studied. Monte 
Carlo simulation is applied to illustrate our proposal. From the result of simulation, 
we find that: (1) The distribution of inclusions in a membrane can be controlled by 
the rate of conformational change of the inclusions. (2) The coupling between inclu- 
sion density and membrane curvature provides the upper limit of the typical size of 
the inclusion clusters. (3) The mobility of inclusions in the membrane depends on 
inclusion-curvature coupling. 

This thesis is organized as follows. In chapter 2 we introduces the model and the 
Hamiltonian of our system. Chapter 3 discusses the simulation algorithm. Chapter 
4 analyzes the result of simulation. Chapter 5 summarizes this thesis. 
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Chapter 2 
The model 



2.1 Introduction 

Since we propose that the formation of heterogeneous structure in biological mem- 
branes can be controlled by the activities (i.e., conformational changes) of the mem- 
brane inclusions, to illustrate our proposal we consider a model in which the mem- 
brane is composed of only one type of lipid molecules and one type of active inclusions 
with two conformational states. These inclusions change their conformational states 
by either external energy drive (ATP, light, etc.), or by ligands. Therefore the mem- 
brane is called an active membrane. Schematics of the conformations of an inclusion 
and its surrounding lipids is shown in Figure ITTl Different conformations of an inclu- 
sion prefer different local membrane curvature, thus the inclusions also induce local 
membrane curvature. We discuss the Hamiltonian that describes our model in the 
following section. 
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Figure 2.1: Schematics of two conformations of an inclusion and its surrounding lipids. 
Different states of an inclusion have different conformation, thus they induce different 
local membrane curvature. 
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2.2 The Hamiltonian 

Our simulation is carried out on a lattice model. In our model, the configuration 
of the system is described by two scalar fields h and n. As shown in Figure 12.21 h 
and n are defined on a square lattice with lattice constant a. The size of a is chosen 
to be the smallest length scale beyond which the continuum elastic description of 
membranes breaks down (a ~ 5nm). [9^ hij is the height of the membrane measured 
from a reference plane (xy-plane) at lattice patch and riij denotes what molecule 
is at patch {i,j) and is defined in the following way, 

, if patch is occupied by lipids, 

1 , if patch is occupied by an inclusion in ground state, (2-1) 

2 , if patch is occupied by an inclusion in excited state. 

The lipids and inclusions are allowed to move in the membrane and the inclusions 

can change their conformations. 

The Hamiltonian of the membrane includes: (1) the bending energy of the mem- 
brane, (2) the surface tension of the membrane, (3) the coupling between inclusion 
density and membrane curvature, and (4) the interaction energy between the inclu- 
sions and lipid molecules. Therefore the discretized form of the hamiltonian of the 
membrane that includes all the above terms can be written as 

.2 



riij = < 



(^i) 

2 

(ij) p=0 

+ ^ Yl (5Z'^™-</'-(^'J>-(^'0) , (2-2) 

((i,j)ik,l)) \m,n J 
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X 

Figure 2.2: The configuration of tlie system is described by h and n. The reference 
plane is discretized into a square lattice with lattice constant a ~ Bnm, and hij is 
the height of the membrane measured from a fiat reference plane at patch A 
white patch is occupied by lipids, and n = at this patch. A gray patch is occupied 
by a ground state inclusion {riij = 1 at this patch). A black patch is occupied by an 
excited state inclusion {riij = 2 at this patch. 
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where k is the membrane bending rigidity constant, and 7 is the surface tension 
constant of the membrane, Cp is curvature-inclusion couphng of type-p component, 
and the last term on the right hand side of Eq. (|2.2|) is the short range interaction 
between the patches, thus Jmn is the interaction energy between one type-m patch 
and a nearest neighboring type-n patch. The short-range interactions between the 
molecules in the membrane could be induced by hydrophobic length mismatch or 
other short range interactions. As shown in Figure 12.31 when the hydrophobic 
length of the inclusions is different from that of a lipid bilayer, the lipid bilayer is 
stretched or compressed from its relaxed thickness. Then free energy of the system is 
increased due to this inclusion-induced thickness variation. Therefore the inclusions 
aggregate together to lower the free energy of the system. Besides the hydrophobic 
length mismatch effect, other short-range interactions between the inclusions also 
exist. means sum over all patches, and J2{{ij){ki)) iiieans sum over all nearest 

neighboring pairs. 0m is defined as follows: 

if riij = 0, then 0o(i, j) = 1, Mhj) = Mhj) = 0, 

iiuij = 1, then j) = 1, Mhj) = Mhj) = 0, (2-3) 

if Uij = 2, then 02(«, j) = 1, Mhj) = Mhj) = 0. 
The discretized bending energy and surface energy are 



and 



a2 



2 



, (2.4) 



2 



(2.5) 
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Figure 2.3: Inclusion-inclusion attraction can be induced by hydrophobic mismatch. 
When the hydrophobic length of the inclusions is different from that of a lipid bilayer, 
the hpid bilayer is stretched or compressed from its relaxed thickness. Then free 
energy of the system is increased due to this inclusion-induced thickness variation. 
Inclusions could aggregate together to lower the free energy of the system. Free energy 
of (b) is lower than (a). 
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Symbol Typical Value 



K ~ 5 X 10^20 . ^ jYY] 

7 ~ 24 X lO^^A^/m [E] 

a ~ 5nm [5] 

/iq ~ Inm 

Cp ~ 1/nm 

Jmn ~ 1 - lO/c^T 

Table 2.1: Typical values of the parameters. 



Typical values of k, 7, a, /iq, Cp, and J^n are listed in Table ITTl For convenience 
the non-dimensionalized hamiltonian is introduced. The unit energy is fcsT, the unite 



length in the xy-plane is a, and the unit length in the z-direction is = ay 
Therefore, the non-dimensionalized hamiltonian is 



lattice 



El 



((i,j)(k,l)) \m,n J 



(2.6) 



{(ij)(fc,0> 

where Gp = hoCp is the non-dimensionalized inclusion-curvature coupling , hij = 



'^^iji Jmn = -j^^Jmni V_L = ^V_L, and 'W\ = 
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Chapter 3 
Simulation method 



Monte Carlo simulation is applied to study the model we introduced in the previous 
chapter. In our simulation, one Monte Carlo step has three parts: (1) The in-plane 
motion of materials in the membrane is simulated by Kawasaki exchange dynam- 
ics. (2) The dynamics of membrane height is simulated by Metropolis algorithm. 
(3) The excitation/relaxation of the inclusions are introduced in the simulation by 
assigning each inclusion certain probability to change its state in each Monte Carlo 
step. When the system has reached the steady state, Hoshen-Kopelman algorithm is 
applied to analyze the size distribution of inclusion clusters in the system. 



3.1 Metropolis algorithm 

Metropolis algorithm provides a set of rules to update the microscopic state of a 
system such that in the long time limit the probability that the system is in a micro- 
scopic state approaches the equilibrium distribution. For example, consider a system 
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described by a single variable x and equilibrium distribution uj{x). Let us assume the 
system to be at state x„ at some instance, we first choose a trial step to change state 
from Xn to Xt- Let 

'J = Uj{Xt)/uj{Xn), (3.1) 

if 7 > 1, then the trial step is accepted and the new state of the system become Xf. 
On the other hand, if 7 < 1, then the trial step has probability 7 to be accepted. 
If the trial step is rejected , then the system remains at state Xn- In the following 
we prove that in the long time limit the probability that the system is in state x 
is given by uj{x). ^3 For an ensemble of independent systems of random initial 
conditions, let be the probability for a system to be in state x at the n-th step, 

and the net probability for a system changing from state x to state y in the next step 



IS 



AP(x -^y)= Pn{x)P{x -^y)- Pn{y)P{y x) 

'P^ix) P{y ^ x) 



Pn{y)P{x ^ y) 



(3.2) 



Where P{x y) is the probability for a system to change its state from x to ?/ in a 
step. When 

Pnjx) ^ Peixl ^ P{y ^ X) 

Pn{y) Pe{y)-P{x^yy ^'^ 
the system is in equilibrium. From equation ()3.2|) . it is clear that if > 

then AP{x ^ y) > 0. On the other hand, when < then AP{x ^ y) < 0. 

Thus, 

Pn{x) Pe{x) 
n^oo P„{y) P^[y) 
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In the Metropolis algorithm, 

P{x ^y) = T{x ^ y)A{x ^ y), (3.5) 

where T{x y) is the probability for a trial step from state x to state y to be chosen, 
A[x ^ y) is the probability for this trial step to be accepted. In the simulation the 
trial steps are chosen such that T{x y) = T{y —>■ x), therefore 

Neix) _ A{y -> x) 



(3.6) 



In Metropolis algorithm, if uj(x) > uj(y), then A(y ^ x) = 1 and A(x ^ y) = ^||y. 
If uj{x) < uj{y), then A{y ^ x) = ^||y and A{x y) = 1. This leads to 

Ne{y) u^iyY ^ ■ ' 

That is, in the long time limit the probability distribution of the system reaches a; (a;). 



3.2 Monte Carlo step 

In our simulation, each Monte Carlo step has 3 parts. The first part is Kawasaki ex- 
change dynamics for the materials in the membrane. We choose one patch randomly, 
and exchange the material on it with a randomly chosen nearest neighbor patch with 
Metropolis rule. From kawasaki exchange dynamics, time scale of one Monte Carlo 
step can be estimated in the following way. The dimension of diffusion constant is 

[D] = Y- (3-8) 

For typical macromolecules diffusing in a membrane, D ~ l/im^/ s PH]- In the absence 
of any interactions, these materials in the membrane move a distance a ~ 5nm per 
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Monte Carlo step, therefore each Monte Carlo step corresponds to a time interval At ~ 
a? /AD ~ 10~^s. In the second part, we choose a new patch at random, and update the 
membrane height at this patch by Ah x ran(— 1, 1) with Metropolis algorithm, where 
Ah is chosen to be ~ OAnm, and ran{—l, 1) is a random number between -1 and 1. 
In the third part, we choose an inclusion randomly, and allow the inclusion to change 
its conformation in the following way: if the inclusion is in the ground (excited) 
state, then it has probability K^n ■ At {K^ff ■ At) to change its conformation to 
excited (ground) state, where {Kgff) is the excitation (relaxation) rate. The first 
and second part are repeated for Niat (total number of patches in the system) times 
per Monte Carlo step. The third part is chosen such that on average each inclusion 
has performed one trial conformational change in each Monte Carlo step. Moreover, 
Metropolis algorithm is applied to the first and the second parts, but in the third 
part the excitation rate and relaxation rate of inclusions are constant. Therefore the 
third part drives the system out of equilibrium, and Kgn, Koff are the important time 
scales which affect the distribution of inclusions. 

The size of the system in the simulations is 128 x 128 to 256 x 256 patches. Total 
number of inclusions in this system is 

A^mc = ^(01 + 02). (3.9) 

For convenience, we also introduce the average inclusion density 0i„c = Ni^c/N. 
The simulations begin with ^(f)incN ground state inclusions and \(t)incN excited state 
inclusions dispersed randomly in a flat discretized membrane. Periodic boundary 
condition is applied in all simulations. After the system has reached steady state, the 
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simulation is performed with up to 10 Monte Carlo steps, and data are taken for 
every 1000 Monte Carlo steps. 

3.3 Statistics of cluster size 

After the system has reached steady state, Hoshen-Kopelman algorithm, a very fast 
and straight method, is applied to analyze the size distribution of inclusion clus- 
ters. |13| Figure 13.11 shows an example of a square lattice containing inclusions (gray 
patches) and lipids (white patches), and the size and number of inclusion clusters are 
to be determined. 

As shown in Figure 13. 2| this is done by going through each row of the lattice 
from left to right, from top to bottom in turn and assigning each gray patch with 
a number n which we call "cluster index" in the following way. The cluster index 
begins from 1. When a gray patch has a nearest-neighbor gray patch that has been 
visited, then the patch is assigned by the same cluster index as this neighbor. On 
the other hand, when a gray patch has no nearest-neighbor gray patch that has been 
visited, then we assign a new number (1 + rimax, where rimax is the maximum n 
that has been assigned up to this patch) to this patch. More than one cluster index 
can be associated with the visited nearest-neighbors of a gray patch. Figure 13.2113.51 
illustrate some examples. The example in Figure 13.21 is when cluster 1 and cluster 
2 are actually the same cluster. In this situation, the Hoshen-Kopelman method 
corrects such "mislabeling" by introducing another set of variable Nn, which is called 
"the label of the n*^^ cluster". If Nn > 0, Nn keeps track of the number of patches 
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Figure 3.2: Labeling patches on a square lattice. This is done by going through each 
row of the lattice from left to right, then from top to bottom in turn and assigning 
each gray patch with a number n (cluster index) in the following way. The cluster 
index begins from 1. If a gray patch has a nearest-neighbor gray patch that has been 
visited, then the patch is assigned by the same cluster index as this neighbor. On the 
other hand, if a gray patch has no nearest-neighbor gray patch that has been visited, 
then we assign a new number {l + rimax, where Umax is the maximum n that has been 
assigned up to this patch) to it. In the third row, more than one cluster indices are 
associated with the nearest neighbors of a gray patch, therefore cluster 1 and cluster 
2 are actually the same cluster. 
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belong to the n cluster. If < 0, denotes that the n cluster actually belongs 
to the cluster labeled by — A^„. 

The way to determine the cluster index in Figure 13.21 is discussed in Figure 13.31 
In Figure ISm the n's and NnS of the nearest neighbors of the new gray patch are 
different, and they belong to different clusters. In Figure 13. 5 1 the n's and NnS of 
the nearest neighbors of the new gray patch are again different, but they actually 
belong to the same cluster. Figure 13.41 and 13.51 explain how Hoshen-Kopelman 
method update both n and A^„. After labeling all patches of the system, we obtain 
the number of clusters and the distribution of cluster size from the final values of 
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Figure 3.3: Hoshen-Kopelman algorithm introduces new variable N^, i.e., "the labels 
of cluster index n". If > 0, is the number of patches belongs to the cluster 
labeled by n. If < 0, denotes that the cluster labeled by n belongs to the 
cluster labeled by —N^- In the right corner of the third row, when there are more 
than one cluster index "n" associated with the nearest neighbors of a new gray patch, 
we check the A^^'s of those neighbors first. In the figure, Nn of the neighboring 
labeled patches are all positive {Ni = 4, N2 = 2). Therefore for the new gray 
patch, n=l (minimum cluster index of the nearest neighbors), and A^i is replaced by 
A'"i + A^2 + 1 = 4 + 2 + 1 = 7, and N2 is replaced by -1, indicating that cluster with 
label n — 2 actually belongs to n = 1 cluster. 
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Figure 3.4: In this case, for the patch labeled by "?" , Nn of its neighboring patches 
are not all positive, A^io = —8, but — —6, Nq — —4, and N4 — 20. In this case, 
the index of the new labeled gray patch is labeled by n = 4, and A^4 is replaced by 
7V4 + iV5 + 1 = 20 + 18 + 1 = 39. A^io and N5 are both replaced by -4. 




Figure 3.5: In this case, Nio = —8, Ng = —6, A'^g = —5, and N5 = 18, thus the 
neighbors of the patch labeled by "?" actually belong to the same cluster. The index 
of this new patch is n = 5, and is replaced by A'^s + 1 = 18 + 1 = 19, and A''io is 
replaced by -5. 
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Chapter 4 

Simulation results and discussion 



4.1 Cluster size distribution 

The simulations begin with ^(pmcN ground state inclusions and |0mc-^ excited state 
inclusions dispersed randomly in a fiat discretized membrane. Periodic boundary 
condition is apphed in all simulations. After the system has reached steady state, the 
simulation is performed with up to 10^ Monte Carlo steps, and data are taken for 
every 1000 Monte Carlo steps. The size of all inclusion clusters can be found by the 
method discussed in Chapter 3. The distribution of inclusion cluster size is defined 
by 



P{uj) = to 



(4.1) 



.Ei!Li(^-^^)J ' 

where P{uj) is the probability that an inclusion is found in a cluster with uj patches, 
and is the total number of clusters with uj patches. Figure shows P{uj) for k = 
5 X 10-20 AT. m, 7 = 24 X IQ-^N/m, 0i„e = 12.5%, K^jfAt = IQ-^ KojKoff = 1/32, 
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Joi = Ji2 = Jo + 1-5, all other J^n = Jo, G2 = 0, Gi = 0, —1, and —2. The maximum 
at = 1 comes from single-inclusion clusters spreading around in the "inclusion-poor" 
domains. Another maximum of P{uj) at greater uj comes from finite-size inclusion 
clusters, this is the typical size of the inclusion clusters. Figure 14.21 shows typical 
snapshots for steady state inclusion distribution and membrane morphology for the 
same parameters as in Figure 14.11 The length scale in the 2;-direction is chosen to be 
the same as x and y directions (the unit length is a) in order to faithfully present the 
morphology of the membrane. One clearly sees that as \Gi\ increases the curvature of 
the membrane in the inclusion-rich domains increases and the typical size of inclusion 
clusters becomes smaller. Furthermore, when Gi 7^ the location of inclusion clusters 
have strong correlation with the regions with high local membrane curvature. This 
is because the system has lower free energy when the inclusions locate in the regions 
with greater membrane curvature. Therefore, when G ^ 0, the membrane prefers to 
form many mountain-like regions with inclusions aggregating on them. 

4.2 Simulation results 

To focus on the mechanisms that we propose in our model clearly, two special cases 
with different interactions between the inclusions and lipids are studied in the simu- 
lations. 

In case 1, J02 = J12 = Jo + J > 0, Joo = Ju = J22 = Joi = Jo, Gi = 0, G2 7^ 0, 
i.e., effectively excited state inclusions tend to attract with each other and induce 
local membrane curvature; ground state inclusions do not attract or repel with other 
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Figure 4.1: P{u!), the probability that an inclusion is found in a cluster with cu- 
patches, for k = 5 x IQ-^^N • m, 7 = 24 x IQ-'^N/m, (pi^c = 12.5%, KffAt = 10'^, 
= -fro///32, Joi = J12 = Jo + 1-5, all other J„„ = Jq, G2 = 0, and Gi = 
(triangle-down), Gi — —1 (square), Gi — —2 (diamond). P{u!) has maximum at 
uj — 1 comes from single-inclusion clusters spreading around in the "inclusion-poor" 
domains. Another maximum of P{u}) at greater M comes from finite-size inclusion 
clusters, this is the typical size of the inclusion clusters and it decreases as |Gi| 
increases due to inclusion-curvature coupling. 



23 




Figure 4.2: Typical snapshots for steady state inclusion distribution and membrane 
morphology for k = 5 x IQ-'^^N-m, 7 = 24 x IQ-^N/m, (pine = 12.5%, K^ffAt = 10-^ 
Kon = Koff/32, Joi = J12 = Jo+ 1.5, all other Jmn = Jo, and G2 = 0. (a) Gi = 0. 
(b) Gi = —1. (c) Gi — —2. As increases, the curvature of the membrane close to 
the inclusion-rich domains increases and the typical size of inclusion clusters becomes 
smaller. The length scale in the 2;-direction is chosen to be the same as x and y 
directions (the unit length is a) in order to faithfully present the morphology of the 
membrane. 
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Figure 4.3: Schematics of inclusion conformation and inclusion-curvature coupling. In 
case 1, the excited state inclusions induce local membrane curvature, and the density 
of ground state inclusions is not coupled to local membrane curvature. In case 2, the 
ground state inclusions induce local membrane curvature, and the density of excited 



state inclusions is not coupled to local membrane curvature. 



molecules, and the density of ground state inclusions is not coupled to local membrane 
curvature. In case 2, Jqi = = Jq + J > 0, Jqo = Ju = J22 = J02 = Jo, G2 = 0, 
Gi 7^ 0, i.e., effectively ground state inclusions tend to attract with each other and 
induce local membrane curvature; excited state inclusions do not attract or repel 
with other molecules, and the density of ground state inclusions is not coupled to 
local membrane curvature. Schematics of the conformations of inclusions are shown 
in Figure 1331 We fix Koff, the relaxation rate of excited inclusions, and study the 
distribution of inclusions for different (j)inc, Kon, Gi, and G2- All simulations are 
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carried out with k = 5 x 10 ^"A^ ■ m, 7 = 24 x 10 ^N/m, a = 5nm, ho = Inm, and 
At ^ 10-^s. 

4.2.1 Short-time in-plane motions of inclusions 

Although the main focus of this thesis is the size distribution of inclusion clusters, 
the effects of conformational change of the inclusions and inclusion-curvature coupling 
to the motion of the inclusions in the membrane are also of great interest. There- 
fore we begin our discussion with short-time in-plan motion of the inclusions. The 
mean square displacement {(P{t)) of an inclusion in short time interval t for case 2 
with (pine = 12.5%, J = 1.5, KoffAt = 10"^, and Kon = Kofj/ 128 is illustrated in 
Figure 14.41 It is clear that {(P) ~ t at sufficiently large t, and the inclusions with 
greater inclusion-curvature coupling move slower in the membrane. Thus we define 
"effective diffusion constant" of an inclusion as the slope of the straight region of the 
(d'^) — t curve. Figure 031 shows the relation between -De// and Kon for case 1 with 
(pine = 12.5%, J = 1.5, and KoffAt = 10'^. It is clear that -Dg// decreases when \G2\ 
increases, but -De// does not depend strongly on Kon- Figure lOl shows the relation 
between -De// and Kon in case 2 with (pine = 12.5%, J = 1.5, and KoffAt = 10~^. It 
is clear that -De// decreases when |Gi| increases, too. Another interesting feature is 
that when Gi = —2, -De// shows strong dependence on Kon- 

Figure ITTI explains the strong dependence of -Dg// on IGI. The free energy is 
lower when a curved inclusion stays in a curved region, and the free energy is higher 
when a curved inclusion stays in a fiat region. Therefore it is less possible for this 
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Figure 4.4: (d^) of an inclusion within short time interval t for case 2 with ^inc — 
12.5%, J = 1.5, Kofft\i = 10-3, = Ko}f/l2S, and = 0. (d^) ~ t at sufficiently 
large t, and the inclusions with stronger inclusion-curvature coupling move slower in 
the membrane. 
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Figure 4.5: D^ff in case 1 with ^i^c = 12.5%, J = 1.5, KoffAt = 10"^, and ^2 = 
(triangle-down), —0.25 (triangle-up) , —0.5 (circle), —1 (square), and —2 (diamond). 
We can find that D^ff decreases as IG2I increases, but D^ff does not depend strongly 
on Kon- At large 16*21, the tending is that Dgff first decreases as increases, then 
the dependence becomes less significant at large Kgn- 
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Figure 4.6: De// in case 2 with 0j„c = 12.5%, J = 1.5, KoffAt = 10"^, and d = 
(triangle-down), —0.25 (triangle-up) , —0.5 (circle), —1 (square), and —2 (diamond). 
It is clear that -De// decreases when IG2I increases. Another interesting feature is that 
when Gi — —2, Deff shows strong dependence on Kgn- 
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XiOW energy 



Figure 4.7: The strong dependence of De// on |G| is explained in this figure. The 
free energy is lower when a curved inclusion stays in a curved region, and the free 
energy is higher when a curved inclusion stays in a flat region. Therefore it is less 
possible for this inclusion to move from a curved region to a flat region. When the 
inclusion-curvature coupling is stronger, this "caging" effect is more significant. 
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inclusion to move from a curved region to a flat region. When the inclusion-curvature 
coupling is stronger, this "caging" effect is more significant. This explains why large 
inclusion-curvature coupling results in small -De//. 

Moreover, in case 2, when Kon increases an inclusion spends less time in the 
ground state. Thus -De// strongly depends on the Kon when \Gi\ is sufficiently large. 
On the other hand, in case 1 Kofj is chosen to be such that Kon > Kofj for all Kon- 
Thus when Kg^ increases, an inclusion spends most of the time in the excited state 
if Kon » Koff. Thus at large 16*215 the tendency is that -De// first decreases as Kon 
increases, then the dependence becomes less significant at large Kon- 

4.2.2 Case 1 

In case 1, when Kof/At > 10^'^ (i.e., Koff > 10^s~^, characteristic time scale of 
a cycle of a motor protein), no inclusion clusters are observed in simulations for 
Kon ^ 10^s~^. This indicates that in a system with typical inclusion excitation 
and relaxation rates, activities may completely suppress the formation of inclusion 
clusters. This is because when Koff is sufficiently large, the lifetime of excited state 
inclusions is sufficiently short such that the excited-state inclusions does not have 
time to aggregate. For Koff At = 10^"^ (i.e., Koff ~ 10s^\ this is more likely for the 
case when the relaxation of an inclusion is induced by ligands) , inclusion clusters are 
observed for a wide range of Kon- 

Figure 14.81 shows the snapshots of the distribution of inclusions in the membrane 
for case 1 with (pine = 12.5%, J = 1.5, and Koff At = 10~^ on a 128 x 128 square 
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lattice. The typical size of inclusion clusters increases as Kon increases or IG2I de- 
creases. Moreover, the dependence of the typical size of inclusion clusters on the 
value of Kon is more evident for small 16*21. Figure shows the relation between 
y/M and Kon^t for KoffAt = 10"^ on a 128 x 128 square lattice. In Figure lOl (a) 
(pine = 12.5%, J = 1.5; in Figure lOl (b) 0j„c = 25%, J = 1. The only exception is 
the empty symbol with the greatest value of v^M in Figure I^Hl (b), which is taken 
from simulations performed on a 256 x 256 square lattice. The filled symbols denote 
the situation when Kon At = the maximum excitation rate can be simulated in our 
simulations. 

Both Figure 14.91 (a) and (b) indicate that when G2 is small, \/M shows weak 
dependence on KonAt for a wide range of i^o^At, but increases fast when KonAt > 

0. 01. This can be understood by comparing the diffusion length of a ground state 
inclusion within its lifetime and the typical size of inclusion clusters. Because the 
diffusion distance during the lifetime of a ground state inclusion is ~ \J AD^ffKon^^- 
Therefore when ^J AD^ffKon^^ ^ a/M, more and more ground state inclusions can 
not escape from an inclusion cluster as Kon increases. This "positive feedback" of 
Kon dependence of v^M is responsible for the strong Kon dependence of \/M at large 
Kon- From the discussion of the previous subsection, -De// ~ O.Sfim'^ / sec ~ 0.8a^/At, 
then the critical length scale is ■\J ADeffKon~^ ~ Kon~^^'^- Therefore, a/M increases 
fast at KonAt > 0.01 {VM ^ 10) in case 1. 

Another important feature of Figure |331 is that for the filled symbols with G2 = 0, 
in Figure EM (a) \/N~'c ~ ^ 45; in Figure |01 (b) V^mc = 64 and a/M ^ 60, 

1. e., almost all inclusions aggregate in the same inclusion cluster . 
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Figure 4.8: The distribution of inclusions in the membrane for case 1 with 0j„c — 
12.5%, J = 1.5, and KgffAt = 10~'^. The typical length scale of inclusion clusters 
increases as Kgn increases or 1^21 decreases. Moreover, the typical size of inclusion 
clusters increases faster when IG2I is smaller. 
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(a): (pine = 12.5%, J = 1.5; (b): (pine = 25%, J = 1. The filled symbols denote 
the situation when KgnAt — 1, the maximum excitation rate can be simulated in 
the simulations, the error bar in this situation is smaller than the symbol. When 
the inclusion-curvature coupling is small, y/M shows weak dependence on KgnAt for 
a wide range of KonAt, but increases fast when KgnAt > 0.01. Furthermore, the 
coupling between inclusion density and membrane curvature provides the upper limit 
of the typical size of the inclusion clusters, and suppresses the formation of large 
inclusion clusters. 
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Moreover, when 6*2 7^ 0, \fM decreases as IG2I increases. This indicates that even 
when all inclusions are all in the excited state, the coupling between inclusion density 
and membrane curvature provides the upper limit of the typical size of the inclusion 
clusters, and suppresses the formation of large inclusion clusters. This is because the 
system has lower free energy when the inclusions locate in the regions with greater 
membrane curvature. Therefore, the membrane prefers to form many mountain-like 
regions with inclusions aggregating on them, when G 7^ 0. 

4.2.3 Case 2 

In case 2, inclusion clusters have been observed for KoffAt = 10~^ — 10~^, and a wide 
range of Kgn- AH simulations are performed on a 128 x 128 square lattice in case 2. 

Figure 101 shows the snapshots of the distribution of inclusions in the membrane 
for case 2 with 0j„c = 12.5%, J = 1.5, and KojfAt = 10^^. The typical size of inclu- 
sion clusters decreases as Kon increases or increases. Moreover, the dependence 
of the typical size of inclusion clusters on the value of Kon is more evident for small 
IGil. Figure E.lll shows the relation between a/M and KgnAt. In Figure E.lll (a) 
^inc = 12.5%, J = 1.5, and K^jfAt = lO'^, in Figure EH] (b) (/)i„c = 25%, J = 1, 
and KoffAt = IQ-^, in Figure EI] (c) 0i„e = 12.5%, J = 1.5, and KfjAt = lO^^. 
The filled symbols denote the situation when all inclusions are in the ground state, 
i.e., Kon = 0. 

When the density of the inclusions is not coupled to local membrane curvature, 
it agrees with a/M ~ Kon~^ relation pretty well for a wide range of Kon- This is 
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Figure 4.10: The distribution of inclusions in the membrane for case 2 with 0j„c = 
12.5%, KoffAt = 10~^, and J = 1.5. The typical length scale of inclusion cluster 
increases as Kon decreases or IGil decreases. Moreover, the typical size of inclusion 
clusters increases faster when IGil is smaller. 
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Figure 4.11: The relation between vM and Kgn in case 1 for (a): 0i„c = 12.5%, 
J = 1.5, and K^ffAt = IQ-^; (b): (pine = 25%, J = 1, and K^ffAt = IQ-^; (c): 
(pine — 12.5%, J = 1.5, KoffAt — 10~^. The filled symbols denote the situation when 
all inclusions are in the ground state, and Kgn = 0, and the error bar in this situation 
is smaller than the symbol. Gi = case agrees with -\/M ~ Kon~^ relation pretty 
well for a wide range of i^on- The slope of the dotted fine is — |. 
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because when there is no budding in the membrane, the growth of inclusion clusters 
in the absence of inclusion activities and inclusion-curvature coupling corresponds to 
a two-dimensional phase separation dynamics, i.e., a/M ~ t^/^. |E1 E| This growth 
eventually saturates due to the active transitions of the inclusions. Since the lifetime 
of ground state inclusions is Kon~^ thus the typical length scales of inclusion clusters 
in the steady state should obey \/M ~ Kon~^^'^- |B] 

The same as case 1, for the filled symbol with Gi = 0, almost all inclusions 
aggregate in the same inclusion cluster (in Figure 14.111 (a) a/A^ a/M «i 45; in 
Figure E. Ill (b) y/N]^ = 64 and v^M ^ 60). Moreover, when Gi ^ 0, v^M decreases 
as |Gi| increases. Therefore, the coupling between inclusion density and membrane 
curvature provides the upper limit of the typical size of the inclusion clusters, and 
suppresses the formation of large inclusion clusters as in case 1. 

Notice that Figure E.lll (a) and (b) have different 0jric and different Jmn, but at 
the same the typical length scales of inclusion clusters are almost in the same 
range. Furthermore, K^jf in (c) is tenfold of -ft'o// in (a), but at the same Kon the 
typical length scales of inclusion clusters are also almost in the same range. This 
means that Kon is the key factor controls the size of inclusion clusters. 

In case 2 the lifetime of the state that prefers to aggregate is changed when K^n 
is changed, but in case 1 the lifetime of the state that prefers to leave the inclusion 
clusters is changed when Kon is changed. Therefore, case 1 and case 2 are different, 
and this is the reason why the Kon dependence in these two cases are so different. 
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Chapter 5 



Summary 



Our Monte Carlo simulations on a simple toy model have shown that inclusion ac- 
tivities and inclusion-curvature coupling both could contribute to the formation of 
finite-size inclusion clusters in biological membranes. Since the inclusions have two 
internal states, and the conformational changes of the inclusions are induced by ex- 
ternal stimuli. Thus, we treat a biomembrane as an active, nonequilibrium system. 
Two special cases with different interactions between the inclusions and lipids are 
studied in the simulation. In case 1, excited state inclusions attract with each other 
and induce local membrane curvature; ground state inclusions do not attract or repel 
with other molecules, and the density of ground state inclusions is not coupled to lo- 
cal membrane curvature. In case 2, ground state inclusions tend to attract with each 
other and induce local membrane curvature; excited state inclusions do not attract or 
repel with other molecules, and the density of excited state inclusions is not coupled 
to local membrane curvature. Our main results are: 

(1) The effective diffusion constant of the inclusions decreases as the inclusion- 
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curvature coupling and the lifetime of the state with inclusion-curvature coupling 
increases. 

(2) In case 1, we find that when Kojf is sufficiently large, the excited-state inclu- 
sions does not have time to aggregate, and there is no inclusion cluster in the system 
for a wide range of Kon- 

(3) In case 1, when inclusion-clusters do form, typical size of clusters shows weak 
dependence on the excitation rate of the inclusions for a wide range of Kon, but 
increases fast when Kon is greater than some critical value. Because the diffusion dis- 
tance of a ground state inclusion in the inclusion cluster ~ AD^ffKon'^ ~ Kon'^^"^- 
When a/M > Kon'^^"^ the inclusions become hard to escape from an inclusion cluster. 

(4) In case 2, we find that when the density of the inclusions is not coupled to 
local membrane curvature, it agrees with a/M ~ Kon~^ relation pretty well for a 
wide range of Kon- This is because when there is no budding in the membrane, 
the growth of inclusion clusters in the absence of inclusion activities and inclusion- 
curvature coupling corresponds to a two-dimensional phase separation dynamics, i.e., 
a/M ~ t^/^. This growth eventually saturates at time scale ~ Kon~^ due to the active 
transitions of the inclusions. 

(5) Moreover, in both case 1 and case 2, we find that the coupling between inclusion 
density and membrane curvature provides the upper limit of the typical size of the 
inclusion clusters, and suppresses the formation of large inclusion clusters. 

Our simple model has neglected many complications that occur in real biological 
membrane. For example: (i) We ignored the effect of cytoskeleton on the motion of the 
inclusions. ^H] (h) The effect of solution flow is not incorporated in our model. JH] (hi) 
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Our model is not appropriate to describe budding in biomembrane. [201 New models 
that includes the aforementioned mechanisms will be our future research direction. 
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APPENDICES 



Appendix A 



Non-dimensionalization of the 



Hamiltonian 



The Hamiltonian of the membrane can be written as 

(id) 

2 

- J] J]a2fi;Cp0p(^,j) (Vi/i).. 

+ 1 Yl (l]'^mn0m(i,J>n(A;,/)) . (A.l) 

First, we choose the unit of energy to be ksT. Then we choose the unite length of 
local vertical height of membrane to be ho ~ Thus the first term of Hiatuce at 

h ho, V_L^ ~ ^ is on the order of /c^T. The discretized Laplacian and gradient of 
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h in two dimension are 



2 



d /'K+q,y - K-^,y^ ^ d fhx,y+^ - hx,y-^ 



dx \ a J dy 

~2 (^^+<^!3/ ~^ hx—a,y ~l~ hx^y-\-a ~l~ hx^y—a '^hx^y) 



0^ I ksT / 



K 



and 



dx dy"^ 



1 

2a 



{hx+a,y hx—a,y^'i ~l~ ihx,y+a hx^y—a)j 

a I ksT 



2a V K 



knT: 



-V±h. 



K 



Thus, total Hamiltonian becomes 



kBTHiattice — 



a 



ij K 



K 



+ ^ ( ^BTJmn(f>mihj)(l>n{k, I) J 

{(i,j){k,l)) \m,n J 
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Prom the definition of dimensionless of inclusion-curvature coupling Gi = hoCi, we 
can obtain the non-dimensionalized Hamiltonian of the system 



Hlnffi, 



lattice 



(hi) L 



2 

ij K 



^3 



(A.5) 
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